---
title: Degradation statistics
subtitle: baRulho:quantifying habitat-induced degradation of (animal) acoustic signals
author: Marcelo Araya-Salas
date: "`r Sys.Date()`"
toc: true
toc-depth: 3
toc-location: left
number-sections: true
highlight-style: pygments
format:
html:
df-print: kable
code-fold: true
code-tools: true
code-copy: true
embed-resources: true
editor_options:
chunk_output_type: console
---
<!-- this code add line numbers to code blocks -->
<!-- only works when code folding is not used in yaml (code_folding: show) -->
```{=html}
<style>
body
{ counter-reset: source-line 0; }
pre.numberSource code
{ counter-reset: none; }
</style>
```
```{r set root directory, echo = FALSE}
# set working directory as project directory or one directory above,
rootdir <- try(rprojroot::find_rstudio_root_file(), silent = TRUE)
if (is(rootdir, "try-error")) rootdir <- ".."
knitr::opts_knit$set(root.dir = rootdir)
```
```{r add link to github repo, echo = FALSE, results='asis'}
# print link to github repo if any
if (file.exists("./.git/config")){
config <- readLines("./.git/config")
url <- grep("url", config, value = TRUE)
url <- gsub("\\turl = |.git$", "", url)
cat("\nSource code, data and annotation protocol found at [", url, "](", url, ")", sep = "")
}
```
```{r setup style, echo = FALSE, message = FALSE, warning=FALSE}
# options to customize chunk outputs
knitr::opts_chunk$set(
class.source = "numberLines lineAnchors", # for code line numbers
tidy.opts = list(width.cutoff = 65),
tidy = TRUE,
message = FALSE,
warning = FALSE
)
```
<div class="alert alert-info">
# Purposes {.unnumbered .unlisted}
- Summarize degradation metrics on re-recorded signals from playback experiment at Bosque de Tlalpan, Mexico City, 2019
- Run statistical analyses
</div>
# Load packages {.unnumbered .unlisted}
```{r load packages}
source("https://raw.githubusercontent.com/maRce10/sketchy/main/R/load_packages.R")
# install/ load packages
load_packages(packages = c("remotes", "kableExtra", "knitr", "formatR", "rprojroot", "maRce10/warbleR", "ggplot2", "tidyr", "viridis", "corrplot", "brms", "ggdist", "cowplot", "cmdstanr", "maRce10/brmsish", "emmeans", "ggsignif"))
my.viridis <- function(...) viridis(alpha = 0.5, begin = 0.3, end = 0.7, ...)
source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
source("~/Dropbox/R_package_testing/brmsish/R/check_rds_fits.R")
```
<!-- this code add line numbers to code blocks -->
<!-- only works when code folding is not used in yaml (code_folding: show) -->
<style>
body
{ counter-reset: source-line 0; }
pre.numberSource code
{ counter-reset: none; }
</style>
```{r read data degradation, eval = TRUE, echo = FALSE, warning=FALSE, message=FALSE, include=FALSE}
degrad_df <- read.csv("./data/processed/barulho_degradation_metrics.csv", stringsAsFactors = FALSE)
degrad_params <- c(bl.rt = "blur.ratio", sp.bl.rt = "spectral.blur.ratio", env.cr = "envelope.correlation", EA ="excess.attenuation", SPL = "SPL", SNR = "signal.to.noise.ratio", SPCC = "cross.correlation", TSR = "tail.to.signal.ratio", TNR = "tail.to.noise.ratio", sp.cr = "spectrum.correlation", pc1 = "pc1")
master_metadata <- read.csv("./data/raw/consolidated.master.sf_selection_table.csv")
# keep only simulated data (exluding data from a student)
master_metadata <- master_metadata[2:961, ]
# infer original frequency
possb_freqs <- seq(0.5, 10, length.out = 20)
master_metadata$freq <- sapply(1:nrow(master_metadata), function(x)
possb_freqs[which.min(abs((master_metadata$top.freq[x] + master_metadata$bottom.freq[x]) / 2 - possb_freqs))]
)
freq_contour <- freq_ts(master_metadata, length.out = 20, parallel = 22, img = FALSE, path = "./data/raw/master_sound_file", pb = FALSE)
freq_contour$mean.freq <- rowMeans(freq_contour[, 3:ncol(freq_contour)])
freq_contour$orig_sound_file <- master_metadata$orig.sound.file
degrad_df$frequency <- sapply(1:nrow(degrad_df), function(x)
freq_contour$mean.freq[freq_contour$orig_sound_file == degrad_df$sound.id[x]]
)
degrad_df$sim.frequency <- sapply(1:nrow(degrad_df), function(x)
master_metadata$freq[master_metadata$orig.sound.file == degrad_df$sound.id[x]]
)
degrad_df <- separate(degrad_df, col = sound.id, into = c("Duration", "Freq_modulation", "Amp_modulation", "Harmonicity", "delete"), remove = FALSE, sep = "_")
degrad_df$Duration <- ifelse(degrad_df$Duration == 0.1, "short", "long")
degrad_df$`Freq_modulation` <- ifelse(degrad_df$`Freq_modulation` == "BB", "fm", "no_fm")
degrad_df$Harmonicity <- gsub(".wav", "", degrad_df$Harmonicity)
degrad_df$delete <- NULL
degrad_df$treatment.replicates <- paste0("frq:", "_", degrad_df$sim.frequency, "dur:", degrad_df$Duration, "_", degrad_df$Amp_modulation, "_", degrad_df$Freq_modulation)
degrad_df$Freq_modulation <- factor(degrad_df$Freq_modulation, levels = c("no_fm", "fm"))
degrad_df$Amp_modulation[degrad_df$Amp_modulation == "no.am"] <- "no_am"
degrad_df$Amp_modulation <- factor(degrad_df$Amp_modulation, levels = c("no_am", "am"))
degrad_df$duration <- factor(degrad_df$Duration, levels = c("short", "long"))
degrad_df <- degrad_df[degrad_df$Harmonicity != "harm", ]
degrad_df <- degrad_df[degrad_df$distance > 1, ]
degrad_df$frequency <- degrad_df$sim.frequency
degrad_df$duration <- degrad_df$Duration
degrad_df$amplitude.modulation <- degrad_df$Amp_modulation
degrad_df$frequency.modulation <- degrad_df$Freq_modulation
degrad_df$location <- degrad_df$Transect
degrad_df$spectrum.blur.ratio <- degrad_df$spectral.blur.ratio
degrad_df$habitat.structure <- ifelse(degrad_df$Vegetacion == "abierta", "open", "closed")
degrad_df$sim.frequency <- degrad_df$signal.to.noise.ratio.t2 <- degrad_df$SPL <- degrad_df$template <- degrad_df$data.set <- degrad_df$Temperatura <- degrad_df$Harmonicity <- degrad_df$freq <- degrad_df$Humedad <- degrad_df$Vegetacion <- degrad_df$Duration <- degrad_df$orig.sound.file <- degrad_df$Freq_modulation <- degrad_df$Amp_modulation <- degrad_df$spectral.blur.ratio <- degrad_df$Transecto <- NULL
names(degrad_df)
degrad_df <- degrad_df[, c("sound.files", "selec", "start", "end", "bottom.freq", "top.freq", "sound.id", "habitat.structure", "distance", "reference", "frequency", "duration", "treatment.replicates", "location", "frequency.modulation", "amplitude.modulation","blur.ratio", "envelope.correlation", "excess.attenuation", "signal.to.noise.ratio", "cross.correlation", "spectrum.blur.ratio", "spectrum.correlation", "tail.to.noise.ratio", "tail.to.signal.ratio")]
write.csv(degrad_df, "./data/processed/tlalpan_degradation_metrics_v01.csv", row.names = FALSE)
```
```{r add PCA, out.width = "100%"}
degrad_df <- read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")
degrad_params <- c("blur.ratio", "spectrum.blur.ratio", "envelope.correlation", "excess.attenuation", "signal.to.noise.ratio", "cross.correlation", "tail.to.signal.ratio", "tail.to.noise.ratio", "spectrum.correlation")
comp.cases <- complete.cases(degrad_df[,names(degrad_df) %in% degrad_params])
pca <- prcomp(degrad_df[comp.cases, names(degrad_df) %in% degrad_params], scale. = TRUE)
# add to data
degrad_df$pc1 <- NA
degrad_df$pc1[comp.cases] <- pca$x[, 1]
degrad_df$pc1.1m.rate <- degrad_df$distance
# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:4])
pca_rot_stck <- stack(pca_rot)
pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)
ggplot(pca_rot_stck, aes(x = variable, y = rotation, fill = Sign)) +
geom_col() +
coord_flip() +
scale_fill_viridis_d(alpha = 0.7, begin = 0.2, end = 0.8) +
facet_wrap(~ ind) +
theme_classic()
```
# Correlation among metrics
Raw metrics:
```{r correlation among raw metrics, out.width = "120%"}
cormat <- cor(degrad_df[, degrad_params], use = "pairwise.complete.obs")
rownames(cormat) <- colnames(cormat) <- degrad_params
cols_corr <- colorRampPalette(c(viridis(3, direction = 1, begin = 0.2, end = 0.5), "#BEBEBE1A", "white", "#BEBEBE1A", viridis(3, direction = 1, begin = 0.7, end = 0.9)))(30)
cp <- corrplot.mixed(cormat, tl.cex = 0.7,
upper.col = cols_corr,
lower.col = cols_corr,
order = "hclust",
lower = "number",
upper = "ellipse",
tl.col = "black")
# sort parameters as in clusters for cross correlation
degrad_params <- degrad_params[match(rownames(cp$corr), degrad_params)]
```
# Data description
- `r length(unique(degrad_df$frequency))` frequencies
- `r length(unique(degrad_df$location))` locations
- `r nrow(degrad_df)` test sounds
- `r length(unique(degrad_df$treatment.replicates))` treatment combinations
Sample sizes per location, transect and signal type:
```{r, results='asis'}
agg <- aggregate(cbind(sound.id, treatment.replicates) ~ location + habitat.structure + distance, degrad_df, function(x) length(unique(x)))
agg$replicates <- agg$sound.id / agg$treatment.replicates
df1 <- knitr::kable(agg, row.names = FALSE,
escape = FALSE, format = "html")
kable_styling(df1, bootstrap_options = c("striped",
"hover", "condensed", "responsive"), full_width = FALSE,
font_size = 15)
```
# Statistical analysis
```{r, eval =FALSE}
iter <- 20000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"), prior(cauchy(0, 4), class = "sd"))
# set frequency to mean-centered and divide by twice the standard deviation
degrad_df$raw_frequency <- degrad_df$frequency
degrad_df$frequency <- (degrad_df$frequency - mean(degrad_df$frequency)) / (2 * sd(degrad_df$frequency))
# set base level for factors
degrad_df$habitat.structure <- factor(degrad_df$habitat.structure, levels = c("open", "closed"))
degrad_df$frequency.modulation <- factor(degrad_df$frequency.modulation, levels = c("no_fm", "fm"))
degrad_df$amplitude.modulation <- factor(degrad_df$amplitude.modulation, levels = c("no_am", "am"))
degrad_df$duration <- factor(degrad_df$duration, levels = c("short", "long"))
degrad_df$location <- as.factor(degrad_df$location)
degrad_df$distance_f <- paste0(degrad_df$distance, "m")
degrad_df$distance_f <- factor(degrad_df$distance_f, levels = c("10m", "30m", "65m", "100m"), ordered = TRUE)
set.seed(123)
cmdstanr::set_cmdstan_path("~/Documentos/cmdstan/")
# to run within-chain parallelization
mod_pc1 <-
brm(
formula = pc1 ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates),
data = degrad_df,
prior = priors,
iter = iter,
chains = chains,
cores = chains,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "./data/processed/regression_model_pc1.RDS",
file_refit = "always",
backend = "cmdstanr",
threads = threading(10)
)
mod_blurratio <-
brm(
formula = blur.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates),
data = degrad_df,
prior = priors,
iter = iter,
chains = chains,
cores = chains,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "./data/processed/regression_model_blur_ratio.RDS",
file_refit = "always",
backend = "cmdstanr",
threads = threading(10)
)
mod_ea <-
brm(
formula = excess.attenuation ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates),
data = degrad_df,
prior = priors,
iter = iter,
chains = chains,
cores = chains,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "./data/processed/regression_model_excess_attenuation.RDS",
file_refit = "always",
backend = "cmdstanr",
threads = threading(10)
)
mod_tsr <-
brm(
formula = tail.to.signal.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates),
data = degrad_df,
prior = priors,
iter = iter,
chains = chains,
cores = chains,
control = list(adapt_delta = 0.99, max_treedepth = 15),
file = "./data/processed/regression_model_tail_to_signal_ratio.RDS",
file_refit = "always",
backend = "cmdstanr",
threads = threading(10)
)
```
## PC1 degradation
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
effects <- c("habitat structure", "frequency modulation", "amplitude modulation", "frequency", "duration", "frequency:habitat structure", "habitat structure:duration", "habitat structure:amplitude modulation", "habitat structure:frequency modulation")
mod <- readRDS("./data/processed/regression_model_pc1.RDS")
extended_summary(mod, n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects, print.name = FALSE, trace = TRUE, return = FALSE)
```
### Contrasts
#### Frequency modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ frequency.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
#### Amplitude modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ amplitude.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
### Conditional plots
```{r, eval = TRUE, warning = FALSE, out.width="100%"}
bs <- 10
ts <- 4
ylm <- c(-5, 1.5)
diff <- 0.1
sd_freq <- sd(degrad_df$frequency)
mean_freq <- round(mean(degrad_df$frequency))
freq_labs <- c(paste(round(mean_freq + sd_freq, 1), "(mean+SD)", sep = "\n"), paste(round(mean_freq - sd_freq, 1), "(mean-SD)", sep = "\n"))
ggf <- plot(conditional_effects(mod, "habitat.structure:frequency", int_conditions = list(frequency = c( mean_freq - sd_freq, mean_freq + sd_freq))), plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = freq_labs) +
theme_classic(base_size = bs) +
labs(color = "Frequency", y = "Degradation (PC1)", x = "Habitat") +
scale_fill_discrete(guide = 'none') +
ylim(ylm)
# add contrasts
for (i in 1:2) ggf <- ggf + geom_signif(y_position = c(1),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
# add larger bracket
ggf <- ggf + geom_signif(y_position = c(1.5),
xmin = c(1), xmax = c(2),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Tonal", "FM")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylm)
# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(0),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggfm <-
ggfm + geom_signif(y_position = c(1),
xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggam <-plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Flat", "AM")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "Degradation (PC1)", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylm)
# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(0),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggam <- ggam + geom_signif(y_position = c(1),
xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Short", "Long")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylm)
# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(0),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggdur <-
ggdur + geom_signif(y_position = c(1),
xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4, color = "gray40")
plot_grid(ggf, ggfm, ggam, ggdur)
```
## Blur ratio
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
mod <- readRDS("./data/processed/regression_model_blur_ratio.RDS")
extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects, print.name = FALSE)
```
### Contrasts
#### Frequency modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ frequency.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
#### Amplitude modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ amplitude.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
### Conditional plots
```{r, eval = TRUE, warning = FALSE, out.width="100%"}
bs <- 10
ts <- 4
ylim <- c(-0.1, 0.3)
ggf <- plot(conditional_effects(mod, "habitat.structure:frequency", int_conditions = list(frequency = c( mean_freq - sd_freq, mean_freq + sd_freq))), plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = freq_labs) +
theme_classic(base_size = bs) +
labs(color = "Frequency", y = "Blur ratio", x = "Habitat") +
scale_fill_discrete(guide = 'none') +
ylim(ylim)
# add contrasts
for (i in 1:2) ggf <- ggf + geom_signif(y_position = c(0.2),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
# add larger bracket
ggf <- ggf + geom_signif(y_position = c(0.23),
xmin = c(1), xmax = c(2),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Tonal", "FM")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
diff <- 0.1
# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(0.25),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggfm <-
ggfm + geom_signif(y_position = c(0.3),
xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggam <-plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Flat", "AM")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "Blur ratio", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(0.2),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggam <- ggam + geom_signif(y_position = c(0.3),
xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Short", "Long")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(0.2),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggdur <-
ggdur + geom_signif(y_position = c(0.3),
xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4, color = "gray40")
plot_grid(ggf, ggfm, ggam, ggdur)
```
## Excess attenuation
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
mod <- readRDS("./data/processed/regression_model_excess_attenuation.RDS")
extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects, print.name = FALSE)
```
### Contrasts
#### Frequency modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ frequency.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
#### Amplitude modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ amplitude.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
### Conditional plots
```{r, eval = TRUE, warning = FALSE, out.width="100%"}
bs <- 10
ts <- 4
ylim <- c(20, 90)
ggf <- plot(conditional_effects(mod, "habitat.structure:frequency", int_conditions = list(frequency = c( mean_freq - sd_freq, mean_freq + sd_freq))), plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = freq_labs) +
theme_classic(base_size = bs) +
labs(color = "Frequency", y = "Excess attenuation", x = "Habitat") +
scale_fill_discrete(guide = 'none') +
ylim(ylim)
# add contrasts
for (i in 1:2) ggf <- ggf + geom_signif(y_position = c(85),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
# add larger bracket
ggf <- ggf + geom_signif(y_position = c(90),
xmin = c(1), xmax = c(2),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Tonal", "FM")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
diff <- 0.1
# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(65),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggfm <-
ggfm + geom_signif(y_position = c(72),
xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggam <-plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Flat", "AM")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "Excess attenuation", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(65),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggam <- ggam + geom_signif(y_position = c(72),
xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Short", "Long")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(65),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggdur <-
ggdur + geom_signif(y_position = c(72),
xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4, color = "gray40")
plot_grid(ggf, ggfm, ggam, ggdur)
```
## Tail-to-signal ratio
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
mod <- readRDS("./data/processed/regression_model_tail_to_signal_ratio.RDS")
extended_summary(fit = mod, n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects, print.name = FALSE)
# mod <- readRDS("./data/processed/regression_model_int.RDS")
```
### Contrasts
#### Frequency modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ frequency.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
#### Amplitude modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ amplitude.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
### Conditional plots
```{r, eval = TRUE, warning = FALSE, out.width="100%"}
bs <- 10
ts <- 4
ylim <- c(-27, -5)
ggf <-
plot(conditional_effects(mod, "habitat.structure:frequency", int_conditions = list(frequency = c( mean_freq - sd_freq, mean_freq + sd_freq))), plot = FALSE)[[1]] + scale_color_viridis_d(end = 0.8, labels = freq_labs) +
theme_classic(base_size = bs) +
labs(color = "Frequency", y = "Tail-to-signal ratio", x = "Habitat") +
scale_fill_discrete(guide = 'none') +
ylim(ylim)
# add contrasts
for (i in 1:2) ggf <- ggf + geom_signif(y_position = c(-9),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
# add larger bracket
ggf <- ggf + geom_signif(y_position = c(-8),
xmin = c(1), xmax = c(2),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggfm <- plot(conditional_effects(mod, "habitat.structure:frequency.modulation"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Tonal", "FM")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
diff <- 0.1
# add contrasts
for (i in 1:2) ggfm <- ggfm + geom_signif(y_position = c(-12),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggfm <-
ggfm + geom_signif(y_position = c(-9),
xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggam <-plot(conditional_effects(mod, "habitat.structure:amplitude.modulation"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Flat", "AM")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "Tail-to-signal ratio", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggam <- ggam + geom_signif(y_position = c(-12),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggam <- ggam + geom_signif(y_position = c(-9),
xmin = 1, xmax = 2,
annotation = c("*"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggdur <- plot(conditional_effects(mod, "habitat.structure:duration"), plot = FALSE)[[1]] +
scale_color_viridis_d(end = 0.8, labels = c("Short", "Long")) +
theme_classic(base_size = bs) +
labs(x = "Habitat", y = "", color = ""
) + scale_fill_discrete(guide = 'none') + theme(legend.position = c(0.93, 0.5)) + ylim(ylim)
# add contrasts
for (i in 1:2) ggdur <- ggdur + geom_signif(y_position = c(-12),
xmin = c(i - diff), xmax = c(i + diff),
annotation = c("ns"), tip_length = 0.02, textsize = ts, size = 0.4, color = "gray40")
ggdur <-
ggdur + geom_signif(y_position = c(-9),
xmin = 1, xmax = 2,
annotation = c("ns"), tip_length = 0.02, textsize = 4, size = 0.4, color = "gray40")
plot_grid(ggf, ggfm, ggam, ggdur)
```
### Contrasts
#### Frequency modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ frequency.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
#### Amplitude modulation * habitat structure
```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"}
cntrs <- as.data.frame(emmeans(mod, pairwise ~ amplitude.modulation * habitat.structure)$contrasts)
names(cntrs)[3:4] <- c("l-95% CI", "u-95% CI")
coef_table <- html_format_coef_table(cntrs, highlight = TRUE, fill = viridis(10)[7])
print(coef_table)
```
## Combined results
```{r, out.width="100%"}
pc1 <- extended_summary(read.file = "./data/processed/regression_model_pc1.RDS", n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)
gg_pc1 <- pc1$plot + labs(x = "Degradation (PC1)")
br <- extended_summary(read.file = "./data/processed/regression_model_blur_ratio.RDS", n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)
gg_br <- br$plot + labs(x = "Blur ratio", y = "") + theme(axis.text.y = element_blank())
ea <- extended_summary(read.file = "./data/processed/regression_model_excess_attenuation.RDS", n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)
gg_ea <- ea$plot + labs(x = "Excess attenuation", y = "") + theme(axis.text.y = element_blank())
tsr <- extended_summary(read.file = "./data/processed/regression_model_tail_to_signal_ratio.RDS", n.posterior = 1000, fill = viridis(10)[7], trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects, print.name = FALSE, trace = FALSE, return = TRUE)
gg_tsr <- tsr$plot + labs(x = "Tail-to-signal ratio", y = "") + theme(axis.text.y = element_blank())
plot_grid(gg_pc1, gg_br, gg_ea, gg_tsr, nrow = 1, rel_widths = c(6, 2, 2, 2))
```
```{r, out.width="100%"}
estimates <- data.frame(pc1 = pc1$coef_table$Estimate, br = br$coef_table$Estimate, ea = ea$coef_table$Estimate, tsr = tsr$coef_table$Estimate)
names(estimates) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation", "Tail-to-signal ratio")
# make estimates relative to maximum estimate in data
rel_estimates <- data.frame(lapply(estimates, function(x) x/ max(abs(x)) * 0.8))
names(rel_estimates) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation", "Tail-to-signal ratio")
# corrplot(as.matrix(estimates), method = "ellipse", )
signif <- data.frame(pc1 = pc1$coef_table$`l-95% CI` * pc1$coef_table$`u-95% CI` > 0, br = br$coef_table$`l-95% CI` * br$coef_table$`u-95% CI` > 0, ea = ea$coef_table$`l-95% CI` * ea$coef_table$`u-95% CI` > 0, tsr = tsr$coef_table$`l-95% CI` * tsr$coef_table$`u-95% CI` > 0)
names(signif) <- c("Degradation (PC1)", "Blur ratio", "Excess attenuation", "Tail-to-signal ratio")
estimates <- as.data.frame(cbind(rownames(pc1$coef_table), stack(estimates)[, 1], stack(rel_estimates)[, 1], stack(signif)))
names(estimates) <- c("predictor", "est", "relavite.est", "sig", "response")
estimates$signif <- ifelse(estimates$sig, "p < 0.05", "p >= 0.05")
# estimates$relavite.est <- ifelse(estimates$sig, estimates$relavite.est, 0)
estimates$response <- factor(estimates$response, levels = rev(c("Degradation (PC1)", "Blur ratio", "Excess attenuation", "Tail-to-signal ratio")))
estimates$predictor <- factor(estimates$predictor, levels = c("habitat structure", "frequency modulation", "amplitude modulation", "frequency", "duration", "frequency:habitat structure", "habitat structure:duration", "habitat structure:amplitude modulation", "habitat structure:frequency modulation"))
ggplot(estimates, aes(x = predictor, y = response, fill = relavite.est)) +
geom_tile() +
coord_equal() +
scale_fill_gradient2(low = viridis(10)[3], high = viridis(10)[7], guide = "none") +
geom_text(aes(label = round(est, 3), color = signif), size = 4) +
scale_color_manual(values = c("black", "gray")) +
labs(x = "", y = "", color = "P value") +
theme_classic(base_size = 10) +
theme(axis.text.x = element_text(color = "black",
size = 11, angle = 30, vjust = 0.8, hjust = 0.8))
```
## Combined model metadata
```{r, warning=FALSE}
check_rds_fits(path = "./data/processed/", html = TRUE)
```
<!-- light green box -->
<div class="alert alert-success">
# Takeaways
- Habitat structure seems to drive most of the observed degradation
- Frequency modulation and amplitude modulation were the signal structural features that more strongly affected transmission
- Signal features interact with habitat structure in diverse ways, sometimes increasing and sometimes decreasing degradation
- Degradation metrics are robust to variation in signal duration
</div>
# Posterior predictive checks
```{r, warning=FALSE, results='asis'}
model_list <- c(
pc1 = "./data/processed/regression_model_pc1.RDS",
blur_ratio = "./data/processed/regression_model_blur_ratio.RDS",
excess_attenuation = "./data/processed/regression_model_excess_attenuation.RDS",
tail_to_signal_ratio = "./data/processed/regression_model_tail_to_signal_ratio.RDS"
)
ndraws <- 20
for (i in seq_len(length(model_list))) {
cat('\n\n## ', names(model_list)[i], '\n\n')
fit <- readRDS(model_list[[i]])
ppc_dens <-
pp_check(fit,
ndraws = ndraws,
type = 'dens_overlay_grouped',
group = "habitat.structure") # shows dens_overlay plot by default
pp_mean <-
pp_check(
fit,
type = "stat_grouped",
stat = "mean",
group = "habitat.structure",
ndraws = ndraws
)
pp_scat <-
pp_check(fit,
type = "scatter_avg_grouped",
group = "habitat.structure",
ndraws = ndraws)
pp_stat2 <- pp_check(fit, type = "stat_2d", ndraws = ndraws)
pp_loo <- pp_check(fit, type = "loo_pit_qq", ndraws = ndraws)
pp_error <-
pp_check(fit,
type = "error_scatter_avg_grouped",
ndraws = ndraws,
group = "habitat.structure")
plot_l <-
list(ppc_dens, pp_mean, pp_scat, pp_error, pp_stat2, pp_loo)
plot_l <-
lapply(plot_l, function(x)
x + scale_color_viridis_d(
begin = 0.1,
end = 0.8,
alpha = 0.5
) + scale_fill_viridis_d(
begin = 0.1,
end = 0.8,
alpha = 0.5
) + theme_classic())
print(plot_grid(plotlist = plot_l, ncol = 2))
}
```
<!-- '---' adds a gray vertical line -->
---
<!-- add packages used, system details and versions -->
<font size="4">Session information</font>
```{r session info, echo=F}
sessionInfo()
```
